#include<stdio.h>
#include <math.h> 
#include<iostream>
#include <direct.h>
#include<ctime>
using namespace std;

//NACA0012翼型生成
double NACA0012(double x)
{
    double y_airfoil;
    y_airfoil=0.6*(-0.1015*pow(x,4)+0.2843*pow(x,3)-0.3576*pow(x,2)-0.1221*x+0.2969*pow(x,0.5));
    return y_airfoil;
}

//简单迭代，返回数组首地址
double* simpleIteration(double **alpha, double **beta, double **gamma, double **ary, int n)
{
    //i为横坐标，j为纵坐标，n为行数，m为列数
    int i,j,m=n+2,k;
    double alpha1[n][m],beta1[n][m],gamma1[n][m],ary1[n][m],ans[n][m],out[2][n][m];
    for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            alpha1[j][i]=*((double*)alpha + m*j + i);
            beta1[j][i]=*((double*)beta + m*j + i);
            gamma1[j][i]=*((double*)gamma + m*j + i);
            ary1[j][i]=*((double*)ary + m*j + i);
        }
    }

    for(j=1;j<n-1;j++){
        for(i=1;i<m-1;i++){
	    ans[j][i] = ((alpha1[j][i]) * (ary1[j][i+1] + ary1[j][i-1]) + \
		gamma1[j][i] * (ary1[j+1][i] + ary1[j-1][i]) - \
		beta1[j][i] / 2.0 * (ary1[j+1][i+1] + ary1[j-1][i-1] - ary1[j-1][i+1] \
		- ary1[j+1][i-1])) / (2.0 * (alpha1[j][i] + gamma1[j][i]));
        }
    }

	for(k=0;k<2;k++){
		for(j=0;j<n;j++){
        	for(i=0;i<m;i++){
				switch(k){
	    	case 0: out[k][j][i] = ans[j][i];
        		}
			}
    	}
	}
    double *pos;
    pos=&out[0][0][0];
	return pos;
}

//寻找juzhen一定范围内最大值
double Find_max(double **ary, int left, int right, int low, int up, int n)
{
	int i,j,m=n+2,k;
    double array[n][m];
    for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            array[j][i]=*((double*)ary + m*j + i);
        }
    }

	double max=array[left][low];
	for(j=low;j<up;j++)
		for(i=left;i<right;i++)
			if(array[j][i]>max) max=array[j][i];
	return max;
}

//O型网格生成，包括翼型、远场、内部、割缝
class Ogrid
{
	//设置访问权限：公开
public:
		int N;
		int radius,iterlimit;
		double errlimit;
		void init( int N_, int radius_, double errlimit_ ,int iterlimit_);
		double* Gridgenerate();
};
	void Ogrid::init( int N_, int radius_, double errlimit_ ,int iterlimit_)
	{
    	N = N_;  radius = radius_; errlimit=errlimit_; iterlimit=iterlimit_;
	}

	double* Ogrid::Gridgenerate()
	{
		int n = N, m=n+2, r = radius;
        double theta[n+2],x[n][m],y[n][m],out1[2][n][m];
        #define pi 3.1415926535898

		int i,j,k;
        for(i=0;i<m;i++){
        theta[i] = 2 * pi * (i - 1) / (n - 1) ;
        }

		// Airfoil boundary
        for(i=1;i<m-1;i++){
		x[0][i] = 0.5 * (1 + cos(theta[i]));
        if (i<(m-1)/2){
		y[0][i] = - NACA0012(x[0][i]);
        }
        else{
		y[0][i] = NACA0012(x[0][i]);
        }
        }

		// Far field boundary
        for(i=1;i<m-1;i++){
		x[n-1][i] = r * cos(-theta[i]);
		y[n-1][i] = r * sin(-theta[i]);
        }

		// In field boundary
		for(j=1;j<n-1;j++){
            for(i=1;i<m-1;i++){
			x[j][i] = x[0][i] + (x[n-1][i] - x[0][i]) * j / (n - 1);
			y[j][i] = y[0][i] + (y[n-1][i] - y[0][i]) * j / (n - 1);
        }
        }

		// Slitting (割缝)
        for(j=0;j<n;j++){
		x[j][0] = x[j][m-3];
        x[j][m-1] = x[j][2];
		y[j][0] = y[j][m-3];
        y[j][m-1] = y[j][2];
        }

    	double alpha[n][m],beta[n][m],gamma[n][m],ary[n][m],x_ks[n][m],x_ita[n][m],y_ks[n][m],y_ita[n][m];
		
		int times=0;
		double err;
		double x_new[n][m], y_new[n][m];
		double *p1,*p2,*p3,*p4;
		p1=&x[0][0];
		p2=&y[0][0];
		do{
			for(j=0;j<n;j++){
        		for(i=0;i<m;i++){
					x_new[j][i]=*((double*)(p1) + m*j + i);
					y_new[j][i]=*((double*)(p2) + m*j + i);
        			}
				}

			// Calculate alpha, beta and gamma
			for(j=1;j<n-1;j++){
        		for(i=1;i<m-1;i++){
			x_ks[j][i] = (x_new[j][i+1] - x_new[j][i-1])*0.5;
			x_ita[j][i] = (x_new[j+1][i] - x_new[j-1][i])*0.5;
			y_ks[j][i] = (y_new[j][i+1] - y_new[j][i-1])*0.5;
			y_ita[j][i] = (y_new[j+1][i] - y_new[j-1][i])*0.5;
				}
			}

			for(j=1;j<n-1;j++){
        		for(i=1;i<m-1;i++){
			alpha[j][i] = pow(x_ita[j][i],2) + pow(y_ita[j][i],2);
			beta[j][i] = x_ks[j][i] * x_ita[j][i] + y_ks[j][i] * y_ita[j][i];
			gamma[j][i] = pow(x_ks[j][i],2) + pow(y_ks[j][i],2);
				}
			}	

			// Using simple iteration equation
			p3=simpleIteration((double**)alpha, (double**)beta, (double**)gamma, (double**)x_new, n);
			for(j=1;j<n-1;j++){
        		for(i=1;i<m-1;i++){
					x[j][i]=*((double*)(p3) + m*j + i);
        			}
				}

			p4=simpleIteration((double**)alpha, (double**)beta, (double**)gamma, (double**)y_new, n);
			for(j=1;j<n-1;j++){
        		for(i=1;i<m-1;i++){
					y[j][i]=*((double*)(p4) + m*j + i);
        			}
				}

			// Slitting
			for(j=0;j<n;j++){
			x[j][0] = x[j][m-3];
			y[j][m-1] = y[j][2];
			x[j][m-1] = x[j][2];
			y[j][0] = y[j][m-3];
			}

			p1=&x[0][0];
			p2=&y[0][0];

			double errx[n][m],erry[n][m];
			for(j=0;j<n;j++){
        		for(i=0;i<m;i++){
			errx[j][i] = abs(x[j][i] - x_new[j][i]);
			erry[j][i] = abs(y[j][i] - y_new[j][i]);
				}
			}
			double errxx=Find_max((double**)errx, 1, m-1, 1, n-1, n);
			double erryy=Find_max((double**)erry, 1, m-1, 1, n-1, n);
			err=max(errxx,erryy);
			times++;
			if (times>iterlimit){
				break;
			}
		}while(err>errlimit);
		
		double out2[2][n][m];
		int u,l;
		for(u=0;u<n;u++){
        	for(l=0;l<m;l++){
				out2[0][u][l]=x[u][l];
				out2[1][u][l]=y[u][l];
        		}
			}
		double *pos2;
        pos2=&out2[0][0][0];

		// Export grid to a dat file
		std::string prefix = "./Output";
		if (_access(prefix.c_str(), 0) == -1)
			_mkdir(prefix.c_str());   
		FILE* fpout;
		fpout = fopen("./Output/O-grid.dat", "w+");
		fprintf(fpout, "variables=\"x\" \"y\"\nzone T=\"O-Grid\" i=%d j=%d F=POINT\n",n,n);
		for (int j = 0; j < n; j++) {
			for (int i = 1; i < m-1; i++){
				fprintf(fpout, "%15.8f\t%15.8f\n", x[j][i],y[j][i]);
			}
		}
		fclose(fpout);
		return pos2;
	}

class GridSolver
{
	//设置访问权限：公开
public:
		int N;
		int radius,iterlimit;
		double errlimit,Vf;
		double *pos,*poss;
		void init( int N_, int radius_, double errlimit_, int iterlimit_, double Vf_ );
		double* phiSolver(double *pos);
		double flowSolver(double *pos, double *poss);
};
	void GridSolver::init( int N_, int radius_, double errlimit_, int iterlimit_, double Vf_)
	{
    	N = N_;  radius = radius_; errlimit=errlimit_; iterlimit=iterlimit_; Vf=Vf_;
	}

	double* GridSolver::phiSolver(double *pos)
    {
    int i,j,n=N,m=n+2;
    double x[n][m],y[n][m],phi[n][m],phii[n][m],phis[n][m],ak[n][m];
    for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            x[j][i]=*((double*)pos + 0*n*m + m*j + i);
            y[j][i]=*((double*)pos + n*m + m*j + i);
    	    }
        }

	for(j=0;j<n;j++){
		for(i=1;i<m-1;i++){
		phi[j][i] = Vf * x[j][i];
		    }
	    }
		
		double alpha[n-2][m-2],beta[n-2][m-2],gamma[n-2][m-2],x_ks,x_it,y_ks,y_it;
		int times=0;
		double err;
		double x_ksi[n], y_ksi[n],x_ita[n],y_ita[n];
		
		do{
			for(j=0;j<n;j++){
			phi[j][0]= phi[j][m-3];
			phi[j][m-1] = phi[j][2];
			}

			double *p1;
			p1=&phi[0][0];
			for(j=0;j<n;j++){
        		for(i=0;i<m;i++){
					phis[j][i]=*((double*)(p1) + m*j + i);
        			}
				}

			x_ks = (x[0][2] - x[0][0]) / 2.0;
			y_ks = (y[0][2] - y[0][0]) / 2.0;
			x_it = (-x[2][1] + 4.0 * x[1][1] - 3.0 * x[0][1]) / 2.0;
			y_it = (-y[2][1] + 4.0 * y[1][1] - 3.0 * y[0][1]) / 2.0;
			phi[0][1] = (4 * phi[1][1] - phi[2][1] - (y_it - x_it) / (y_ks - x_ks) * (phi[0][2] - phi[0][0])) / 3.0;
			phi[0][-1] = phi[0][1];

			for(i=0;i<n;i++){
			x_ksi[i] = (x[0][i+3] - x[0][i+1]) / 2.0;
			y_ksi[i] = (y[0][i+3] - y[0][i+1]) / 2.0;
			x_ita[i] = (-x[2][i+2] + 4.0 * x[1][i+2] - 3.0 * x[0][i+2]) / 2.0;
			y_ita[i] = (-y[2][i+2] + 4.0 * y[1][i+2] - 3.0 * y[0][i+2]) / 2.0;
			}

			for(i=0;i<m-4;i++){
			phii[0][i+1] = (phi[0][i+3] - phi[0][i+1]) / 2.0;
			ak[0][i+1] = (x_ksi[i] * x_ita[i] + y_ksi[i] * y_ita[i]) / (pow(x_ksi[i],2) + pow(y_ksi[i],2));
			phi[0][i+2] = - (2.0 * phii[0][i+1] * ak[0][i+1] - 4.0 * phi[1][i+2] + phi[2][i+2]) / 3.0;
			}

			for(j=1;j<n-1;j++){
			phi[j][0] = phi[j][m-3];
			phi[j][m-1] = phi[j][2];
			}

			for(j=0;j<n-2;j++){
				for(i=0;i<m-2;i++){
					alpha[j][i] = (pow((x[j+2][i+1] - x[j][i+1]),2) + pow((y[j+2][i+1] - y[j][i+1]),2)) / 4.0;
					beta[j][i] = -(((x[j+2][i+1] - x[j][i+1]) * (x[j+1][i+2] - x[j+1][i])) + \
						((y[j+2][i+1] - y[j][i+1]) * (y[j+1][i+2] - y[j+1][i]))) / 4.0;
					gamma[j][i] = (pow((x[j+1][i+2] - x[j+1][i]),2) + pow((y[j+1][i+2] - y[j+1][i]),2)) / 4.0;
					}
			}

			for(j=0;j<n-2;j++){
				for(i=0;i<m-2;i++){
					phi[j+1][i+1] = 0.5 * (alpha[j][i] * (phi[j+1][i+2] + phi[j+1][i]) + gamma[j][i] * \
					(phi[j+2][i+1] + phi[j][i+1]) + 0.5 * beta[j][i] * (phi[j+2][i+2] + phi[j][i] -\
					phi[j][i+2] - phi[j+2][i])) / (alpha[j][i] + gamma[j][i]);
				}
			}

			double errr[n][m];
			for(j=0;j<n;j++){
        		for(i=0;i<m;i++){
			errr[j][i] = abs(phi[j][i] - phis[j][i]);
				}
			}
			err=Find_max((double**)errr, 1, m-1, 1, n-1, n);
			times++;
			if (times==iterlimit){
				break;
			}
		}while(err>errlimit);
		double *pos3;
        pos3=&phi[0][0];
		return pos3;
	}

	double GridSolver::flowSolver(double *pos,double *poss)
	{
	int i,j,n=N,m=n+2;
    double x[n][m],y[n][m],phi[n][m],u[n][m],v[n][m],Cp[n][m],Jacobi[m-3],Jacobi2[n-2][m-2];
    for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            x[j][i]=*((double*)pos + 0*n*m + m*j + i);
            y[j][i]=*((double*)pos + n*m + m*j + i);
    	}
    }
    for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            phi[j][i]=*((double*)poss + 0*n*m + m*j + i);
    	}
    }

	for(i=0;i<m-3;i++){
		Jacobi[i] = (x[0][i+3] - x[0][i+1]) * (-3.0 * y[0][i+2] + 4.0 * y[1][i+2] - y[2][i+2]) -\
			(y[0][i+3] - y[0][i+1]) * (-3.0 * x[0][i+2] + 4.0 * x[1][i+2] - x[2][i+2]);
	}

	for(i=0;i<m-3;i++){
		u[0][i+2] = ((phi[0][i+3] - phi[0][i+1]) * (-3.0 * y[0][i+2] + 4.0 * y[1][i+2] - y[2][i+2]) -\
			(-3.0 * phi[0][i+2] + 4.0 * phi[1][i+2] - phi[2][i+2]) * (y[0][i+3] - y[0][i+1])) / Jacobi[i];
		v[0][i+2] = ((x[0][i+3] - x[0][i+1]) * (-3.0 * phi[0][i+2] + 4.0 * phi[1][i+2] - phi[2][i+2]) -\
			(-3.0 * x[0][i+2] + 4.0 * x[1][i+2] - x[2][i+2]) * (phi[0][i+3] - phi[0][i+1])) / Jacobi[i];
		Cp[0][i+2] = 1 - (pow(u[0][i+2],2) + pow(v[0][i+2],2)) / pow(Vf,2);
	}
		u[0][1] = 0.0;
		v[0][1] = 0.0;
		Cp[0][1] = 1.0;
		for(i=1;i<m-1;i++){
		u[n-1][i] = Vf;
		v[n-1][i] = 0.0;
		Cp[n-1][i] = 0.0;
		}
	for(j=0;j<n-2;j++){
		for(i=0;i<m-2;i++){
		Jacobi2[j][i] = (x[j+1][i+2] - x[j+1][i]) * (y[j+2][i+1] - y[j][i+1]) - \
			(y[j+1][i+2] - y[j+1][i]) * (x[j+2][i+1] - x[j][i+1]);
		u[j+1][i+1] = ((phi[j+1][i+2] - phi[j+1][i]) * (y[j+2][i+1] - y[j][i+1]) -\
			(y[j+1][i+2] - y[j+1][i]) * (phi[j+2][i+1] - phi[j][i+1])) / Jacobi2[j][i];
		v[j+1][i+1] = ((x[j+1][i+2] - x[j+1][i]) * (phi[j+2][i+1] - phi[j][i+1]) -\
			(phi[j+1][i+2] - phi[j+1][i]) * (x[j+2][i+1] - x[j][i+1])) / Jacobi2[j][i];
		Cp[j+1][i+1] = 1 - (pow(u[j+1][i+1],2) + pow(v[j+1][i+1],2)) / pow(Vf,2);
		}
	}
   
		FILE* fpout_flowfiled;
		fpout_flowfiled = fopen("./Output/Flowfield.dat", "w+");
		fprintf(fpout_flowfiled, "variables=\"x\" \"y\" \"u\" \"v\" \"Cp\"\nzone T=\"O-Grid-Flowfield\" i=%d j=%d F=POINT\n",n,n);
		for (int j = 0; j < n; j++) {
			for (int i = 1; i < m-1; i++){
				fprintf(fpout_flowfiled, "%15.8f\t%15.8f\t%15.8f\t%15.8f\t%15.8f\t\n", x[j][i],y[j][i],u[j][i],v[j][i],Cp[j][i]);
			}
		}
		fclose(fpout_flowfiled);

		FILE* fpout_Cp;
		fpout_Cp = fopen("./Output/Cp_x.dat", "w+");
		fprintf(fpout_Cp, "variables=\"x\" \"Cp\"\nzone T=\"O-Grid-Cp\" i=%d F=POINT\n",n-2);
		for (int i = 1; i < m-3; i++){
			fprintf(fpout_Cp, "%15.8f\t%15.8f\t\n", x[0][i],Cp[0][i]);
		}
		fclose(fpout_Cp);
		return 0;
	}

int main()
{
	clock_t startTime,endTime;
	startTime = clock();//计时
	printf("The program is running...\n");

    //从Input.txt文档输入参数
	FILE *fp = NULL;
    char *file = (char*)"input.txt";//input文档路径
    char line[30];
    if( (0 != access(file,R_OK|F_OK)) || (NULL==(fp=fopen(file,"r"))) )
    {
        printf("open %s failed\n",file);
        return -1;
    }
    while( fgets(line, 256, fp) != NULL ){}
    if(fp!=NULL)
    {
        fclose(fp);
    }
	const char *d = ",";
    char *p;
	int kk=0;
	double output[5];
    p = strtok(line,d);
    while(p)
    {
		char str[10];
		string ing=p;
		int ii;
    for( ii=0;ii<ing.length();ii++)
        str[ii] = ing[ii];
    str[ii] = '\0';
		double dd;
		output[kk]=atof(str);
        p=strtok(NULL,d);
		kk++;
    }

	int n=int(output[0]),r=int(output[1]),iterlimit=int(output[3]);//n为径向节点数，r为半径
	double errlimit=output[2],Vf=output[4];//errlimit为误差限，Vf为来流速度
    int i,j,k,m=n+2;
	double *pos_1,*pos_2,*pos_3,*pos_4;
	Ogrid ogrid;
	ogrid.init(n,r,errlimit,iterlimit);
	pos_1=ogrid.Gridgenerate();//生成网格，返回x,y
	printf("O-Grid Generation  Fineshed    step 1 in 3\n");

	double out[2][n][m];
	for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            out[0][j][i]=*((double*)pos_1 + 0*n*m + m*j + i);
            out[1][j][i]=*((double*)pos_1 + n*m + m*j + i);
    	}
    }
	pos_2=&out[0][0][0];
	GridSolver solver;
	solver.init( n, r, errlimit, iterlimit, Vf);
	pos_3=solver.phiSolver(pos_2);//计算phi

	double phi[n][m];  
	for(j=0;j<n;j++){
        for(i=0;i<m;i++){
            phi[j][i]=*((double*)pos_3 + 0*n*m + m*j + i);
    	}
    }
	pos_4=&phi[0][0];
	printf("Grid Solver        Fineshed    step 2 in 3\n");
	solver.flowSolver(pos_2,pos_4);//计算流场
	printf("Flowfiled Solver   Fineshed    step 3 in 3\n");

	endTime = clock();
	printf("Duration=%7.3fs\n",(double)(endTime - startTime)/1000);
	system("pause");
    return 0;
}
